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Abstract: We study the solution of the diffusion equation for Ultra-High Energy Cosmic Rays in the 
general case of an expanding universe, comparing it with the well known Syrovatsky solution obtained 
in the more restrictive case of a static universe. The formal comparison of the two solutions with all 
parameters being fixed identically reveals an appreciable discrepancy. This discrepancy is less important 
if in both models a different set of best-fit parameters is used. 
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Introduction 

Diffusive propagation of Ultra High Energy Cos- 
mic Rays (UHECR) in extragalactic space has been 
recently studied by [1, 2, 8, 3] using the Syrovatsky 
solution (see [9]) of the diffusion equation. The 
Syrovatsky solution is obtained under the restric- 
tive assumptions of time-independent diffusion co- 
efficient (D = D{E)) and energy losses of par- 
ticles (dE/dt = b(E)). Recently two papers ap- 
peared [5, 6] solving the problem of the generaliza- 
tion of the diffusion equation (and its solution) in 
the case of an expanding universe, i.e. in the case 
of time dependent diffusion coefficient and energy 
losses. In these works an analytic solution of the 
diffusion equation in an expanding universe was 
found, valid in the general case of time-dependent 
diffusion coefficient and energy losses, we will re- 
fer to this solution as the Berezinsky-Gazizov (BG) 
solution [5]. In the present paper we will compare, 
following the approach of [6], the spectra com- 
puted in the generalized case (BG solution) and 
the spectra obtained with the Syrovatsky solution 
as in the above cited papers. The diffusion equation 
for ultra-relativistic particles propagating in an ex- 
panding universe from a single source, as obtained 
in [5], reads 
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where the coordinate x corresponds to the comov- 
ing distance and a(t) is the scaling factor of the 
expanding universe, n — n(t, x, E) is the parti- 
cle number density per unit energy in an expanding 
volume, dE/dt = —b(E, t) describes the total en- 
ergy losses, which include adiabatic H(t)E as well 
as interaction b int (E, t) energy losses. Q S (E, t) is 
the generation function, that gives the number of 
particles generated by a single source at coordinate 
x g per unit energy and unit time. 

According to [5], the spherically-symmetric solu- 
tion of Eq. (1) is 
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with cosmological parameters fi m = 0.27 and 
A = 0.73, 



, z) = [ dz' 
Jo 



dt' 



dz 1 



D(E g ,z') 
a 2 (z') ' 



(4) 



UHECR Diffusion in an Expanding Universe 



exp 



dE g (E, z) 
dE 

dz 1 



dt' 



dz 1 



= (1 + *) 
db int (E g ,z') 



dE n 



(5) 



The generation energy E g — E g {E, z) is the solu- 
tion of the energy-losses equation: 



dEg 
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-[H(t)E g + b int (E g ,t)] 
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with initial condition E g (E, 0) = E. 

In the present paper we will discuss the propa- 
gation of UHE protons in Intergalactic Magnetic 
Fields (IMF) following the approach used by [1, 2], 
in which the IMF is produced by a turbulent mag- 
netized plasma. In this picture the IMF is charac- 
terized by a coherent field B c on scales I > l c , 
where l c is the basic scale of turbulence, and on 
smaller scales I < l c the IMF is determined by its 
(assumed) turbulent spectrum. In our estimates we 
will keep l c ~ 1 Mpc. 

The propagation of UHE protons in IMF is char- 
acterized by two basic scales: an energy scale 
E c that follows from the condition r^(i? c ) = l c , 
with tl Larmor radius of the proton, and the dif- 
fusion length ld(E), that is defined as the dis- 
tance at which a proton is scattered by 1 rad. Us- 
ing ld(E) the diffusion coefficient is defined as 
D(E) = cl d (E)/3. 

We can easily identify two separate regimes in the 
particle propagation in IMF, that follows from the 
comparison of the two scale n, and l c . In the case 
r L (E) > l c (E > E c ) the diffusion length can 
be straightforwardly found from multiple scatter- 
ing as 



l d (E) = 1.2 



E 18 
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Mpc, 
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where E 18 = E/(10 18 eV) and B nG = B/{1 nG). 
At E — E c , Id — l c - In the opposite scenario when 
tl < l c (E < E c ) the diffusion length depends on 
the IMF turbulent spectrum. In this case, following 
[6], we have assumed two different pictures: the 
Kolmogorov spectrum l d (E) = l c (E /Ec) 1 / 3 and 
the Bohm spectrum l d (E) = l c (E/E c ). 

The strongest observational upper limit on the IMF 
in our picture is given by [7] as B c < 10 nG on the 
turbulence scale l c = 10 Mpc. In the calculations 
presented here we assume a typical value of B c in 
the range (0.1 - 1) nG and l c = 1 Mpc. 



In the present paper we will not perform a de- 
tailed discussion of the proton diffusion in the gen- 
eral case of an expanding universe, we will ad- 
dress this issue in a forthcoming paper [4], our 
main goal here is to perform a detailed compar- 
ison of the BG solution with the Syrovatsky solu- 
tion. As already discussed in [6], the difference be- 
tween these two solutions is substantial at energies 
E < 3 x 10 18 eV, where the effect of the universe 
expansion (in particular, of the CMB temperature 
growth with red-shift) is not negligible. The high 
energy tail of the UHECR spectrum is less affected 
by the expansion of the universe, nevertheless it is 
interesting to test the compatibility of the BG and 
Syrovatsky spectra at these energies where a sub- 
stantial agreement of the two is expected. 

Diffusive energy spectra of UHECR 

In the present calculations we used a simplified 
description of the IMF evolution with redshift, 
namely we parametrize the evolution of magnetic 
configuration (l c ,B c ) as 

l c (z) = l c /(l + z), B c (z) = B C (1 + z) 2 - m , 

where the term (1 + z) 2 describes the depletion of 
the magnetic field with time due to the magnetic 
flux conservation and (1 + z)~ m due to MHD am- 
plification of the field. The critical energy E c (z) 
found from r L (E) = l c {z) is given by 



E c (z) = 0.93 x 10 18 (1 + 



B c 
1 nG 



for l c = 1 Mpc. The maximum redshift used in the 
calculations is z max = 4. 

Following [2], we have computed the diffuse flux 
assuming a distribution of sources on a lattice with 
spacing d and an injection spectrum, equal for all 
sources, given by 
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where E is a normalizing energy (we used E = 
1 x 10 18 eV) and q n represents the source luminos- 
ity in protons with energies E > Eq, L p (> Eq). 
The corresponding emissivity Co = qo/d 3 , i.e. the 
energy production rate in particles with E > E 
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Figure 1 : Convergence of the diffusive solution to 
the universal spectrum when the distance between 
sources diminishes from 50 to 10 Mpc shown by 
numbers on the curves. 

per unit comoving volume, will be used to fit the 
observed spectrum by the calculated one. 

In figure 1 we test the BG solution with the help 
of the diffusion theorem [1], which states that the 
diffusive solution converges to the universal spec- 
trum, i.e. the flux computed with rectilinear propa- 
gation for an homogeneous distribution of sources, 
in the limit d — > 0, being d the lattice spacing. Fig- 
ure 1 clearly shows this convergence even in the 
case of a strong magnetic field B c = 100 nG (and 
Kolmogorov diffusion). 

In the case of a small distance between source 
and observer the diffusive approximation is not 
valid. This result follows from a simple argument, 
the diffusive approximation works if the diffusive 
propagation time r 2 /D is larger than the time of 
rectilinear propagation, r/c. This condition, us- 
ing D ~ eld, results in r > Id- At distances 
r < Id the rectilinear and diffusive trajectories in 
IMF differ by a little quantity and rectilinear prop- 
agation is a good approximation as far as spec- 
tra are concerned. The number densities of par- 
ticles Q/Aircr 2 and Q/AnDr, calculated in recti- 
linear and diffusive approximations, respectively, 
are equal at r ~ Id, where Q is the rate of parti- 
cle production. We calculated the number densities 
of protons n(E, r) numerically for both modes of 
propagations with energy losses of protons taken 



Figure 2: Equal parameter comparison of the BG 
(expanding universe) and Syrovatsky (static uni- 
verse) solutions, for j g = 2.7, £o = 2.4 x 
10 45 erg/Mpc 3 yr and d — 30 Mpc. The magnetic 
field configuration assumed is B c — 0.1 nG and 
l c = 1 Mpc with different diffusion regimes as in- 
dicated on the plot. 

into account, and the transition is taken from the 
equality of the two spectra. We know that this 
recipe is somewhat rough and an interpolation be- 
tween the two regimes is required [2]. However, 
this interpolation is somewhat difficult because the 
diffusive regime sets up at distances not less than 
six diffusion lengths Id- At distances Id < r < 6 Id 
some intermediate regime of propagation is valid. 
When studied in numerical simulations (e.g. [10]), 
the calculated number density n(E, r) satisfies the 
particle number conservation Anr 2 nu = Q, where 
u is the streaming velocity, while with a simple in- 
terpolated spectrum this condition is not fulfilled a 
priori. In the present paper we will not address this 
problem, that will be studied in a forthcoming pa- 
per [4], assuming the rough recipe for the transition 
between diffusive and rectilinear regimes depicted 
above. This computation scheme can produce ar- 
tificial features in the spectra, that are useful as a 
mark of the transition between the two regimes. 

The direct comparison of the BG and Syrovatsky 
solutions of the diffusion equations is not possi- 
ble because they are embedded in different cos- 
mological environments. While the BG solution is 
valid for an expanding universe, the Syrovatsky so- 
lution is valid only for a static universe. Using two 
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Figure 3: Best fit comparison of the BG (expand- 
ing universe) and Syrovatsky (static universe) so- 
lutions, for 7 9 = 2.7, C = 2.4 x 10 45 erg/Mpc 3 yr 
and d — 30 Mpc. The magnetic field configura- 
tion assumed is B c = 0.1 nG and l c = 1 Mpc with 
different diffusion regimes as indicated on the plot. 



different cosmological models for these solutions, 
there are two ways of comparison. The first one is 
given by equal values of parameters in both solu- 
tions. In this method for BG solution we use the 
standard cosmological parameters for an expand- 
ing universe Hq, f2 TO , A and maximum red-shift 
Zmax U P to which UHECR sources are still active, 
magnetic field configuration (B c ,l c ), separation d 
and UHECR parameters j g and Cq, determined by 
the best fit of the observed spectrum. For a static 
universe with Syrovatsky solution we use the same 
parameters H , d, {B c , l c ), 7 g and Cq. The max- 
imum red-shift in the BG solution is fixed by the 
age of the universe which equals to to = H^ 1 
in the static universe {z max — 1.5). This for- 
mal method of comparison will be referred to as 
"equal-parameter method". Physically a better jus- 
tified comparison is given by the best fit method, in 
which j g and Co are chosen as the best fit parame- 
ters for both solutions independently. 

The comparison of the two solutions is given in 
Figures 2 and 3 in the case of B c = 0.1 nG 
and l c = 1 Mpc with a source spacing d = 30 
Mpc. From these figures one can see a reason- 
ably good agreement between the Syrovatsky solu- 
tion, embedded in a static universe model, with the 
BG solution for an expanding universe at energies 



E > 1 x 10 18 eV, at smaller energies appears a 
noticeable discrepancy between the two solutions 
that is natural and understandable as discussed in 
the introduction. We conclude stating that, from a 
physical point of view, the second method of com- 
parison is more meaningful and it gives a substan- 
tial agreement of the spectra obtained in the two 
cases. 
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